Capturing correlations in chaotic diffusion by approximation methods 
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We investigate three different metiiods for systematically approximating the diffusion coefficient of a deter- 
ministic random walk on the line which contains dynamical correlations that change irregularly under parameter 
variation. Capturing these correlations by incorporating higher order terms, all schemes converge to the analyt- 
ically exact result. Two of these methods are based on expanding the Taylor-Green-Kubo formula for diffusion, 
whilst the third method approximates Markov partitions and transition matrices by using a slight variation of 
the escape rate theory of chaotic diffusion. We check the practicability of the different methods by working 
them out analytically and numerically for a simple one-dimensional map, study their convergence and critically 
discuss their usefulness in identifying a possible fractal instability of parameter-dependent diffusion, in case of 
dynamics where exact results for the diffusion coefficient are not available. 
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I. INTRODUCTION 

Diffusion is a fundamental macroscopic transport process 
in many-particle systems. It is quantifiable by the diffusion 
coefficient, which describes the linear growth in the mean- 
squared displacement of an ensemble of particles. The source 
of this growth is often considered to be a Brownian or ran- 
dom process of collisions between particles. However, on a 
microscopic scale the equations governing these collisions in 
physical systems are deterministic and typically chaotic. By 
studying diffusion in chaotic dynamical systems we can at- 
tempt to take these deterministic rules into account and un- 
derstand the phenomenon of diffusion from first principles 1 1- 
|4t]. Of particular interest is the study of the diffusion coeffi- 
cient under parameter variation in chaotic dynamical systems 
such as one-dimensional maps IJSH^, area preserving two di- 
mensional maps ll9l- [ll]l and particle billiards II12I - U5I1 . Where 
exact analytical results for chaotic dynamical systems exist 
llla42lll . one finds that the diffusion coefficient is typically 
a complicated fractal function of control parameters. This 
phenomenon can be understood as a topological instability of 
the deterministic diffusive dynamics under parameter varia- 
tion il[iMl. 

So far exact analytical solutions for the diffusion coefficient 
could only be derived for simple cases of low-dimensional dy- 
namics. In higher dimensions even very fundamental prop- 
erties of diffusion coefficients are often unknown, such as 
whether they are smooth or fractal functions of control pa- 
rameters Jj, [Tj, 1221 ■ For example, much effort was spent two 
decades ago studying more complicated systems like a two- 
dimensional family of sawtooth maps i l0ll23ll24tl . However, 
despite a good understanding of the orbit structure 11251 12611 it 
was not possible to conclude whether the diffusion coefficient 
is fractal or not [27]. If one wishes to achieve a microscopic 
understanding of diffusion in more realistic physical systems, 
one therefore has to rely either on numerical simulations or on 
approximation methods. 
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In this paper we compare three different methods for ap- 
proximating parameter dependent diffusion coefficients with 
each other by working them out analytically and numerically 
for a simple one-dimensional map. This model has the big 
advantage that it is very amenable to rigorous anal ysis . Its 
diffusion coefficient has been calculated exactly in flV] and 
was found to be a fractal function of a control parameter. Our 
goal is to assess the individual capabilities and limitations of 
these approximation methods in terms of practicability, phys- 
ical interpretation, convergence towards the exact result, and 
identification of an underlying fractal structure in the diffu- 
sion coefficient. We also address recent criticism by Gilbert 
and Sanders [28], who claimed that one of these methods, as 
originally proposed in |29], is mathematically wrong and un- 
physical. 

In sectionlUwe define the deterministic dynamical system 
that provides our test case, which is a simple piecewise linear 
one-dimensional map. In section HU] the first approximation 
method is introduced, called coiTelated random walk in 112911 . 
which consists of truncating the Taylor-Green-Kubo formula 
for diffusion. This method enables us to analytically build up 
a series of approximations which gives evidence for a fractal 
structure. In previous work this approximation scheme has 
successfully been applied to understand parameter dependent 
diffusion in models that are much more complicated than the 
one considered here L3^ ^ J4, 151 . Motivated by the criti- 
cism of ll28ll . in this paper we provide further insight into 
the functioning of this method by working it out rigorously 
for our specific example. In section |IV] the persistent ran- 
dom walk method for diffusion is studied. This method was 
originally proposed within stochastic theory in the form of a 
persistent random walk 1301 13111 . It consists of approximat- 
ing the Taylor-Green-Kubo formula by including memory in 
a self-consistent, persistent way. Recently this method has 
been worked out for chaotic diffusion in Hamiltonian parti- 
cle billiards 11281 l32l 13311 . Here we apply this scheme to the 
different case of a one-dimensional map, and we obtain a se- 
ries of approximations analytically and then numerically. In 
section [V] we look at a third method, defined by a slight vari- 
ation 16-18] of the escape rate theory of chaotic diffusion 
ill I2I I34I437I1 in that absorbing boundary conditions are re- 
placed by periodic ones. This method thus consists of eval- 



uating the diffusion coefficient in terms of the decay rate of 
the dynamical system towards the equilibrium state, instead 
of using the escape rate. The decay rate is in turn obtained 
by an approximation to the relevant Markov transition matrix. 
By this method we are able to build up a series of approxima- 
tions which, through the functional form of the interpolation 
that we find, gives very strong evidence for fractahty. Basic 
ideas defining this method have been sketched in yl, how- 
ever, this is the first time that it has been fully worked out to 
understand fractal diffusion coefficients. Section[Vl]forms the 
conclusion. 



II. A ONE-DIMENSIONAL MAP EXHIBITING CHAOTIC 
DIFFUSION 

We use the simplest setting possible, where determinis- 
tic diffusion is generated by a parameter dependent one- 
dimensional dynamical system. The equations of motion are 
determined by a map Mh (x) : R ^> M so that 



Xn+l = Mh{x„) 



M]l+\x) x£ 



l,h>0,nen, 



(1) 



with X = xo IdhTI [iTIl . In our case, the map Mh{x) is based 
on the Bernoulli shift or doubling map, combined with a lift 
parameter h, which gives the simple parameter dependent map 
of the interval 



Mh{x) = 



2x + h 0<x<i 
2x-l-h i<2:<l 



(2) 



This map exhibits 'escape', i.e., points leave the unit interval 
under iteration. It is copied and lifted over the real line by 



Mh{x + z) = Mh{x) + z, z eZ 



(3) 



in order to obtain a map from the real line to itself, see 
Fig- Ea)- The symmetry in this system ensures that there 
is no mean drift 1 19]. Note that the invariant density of the 
map Eq. (|2]i modulo 1 remains by construction simply uni- 
form throughout the whole parameter range. This is in con- 
trast to the related piecewise linear maps studied in ll3l [l6l - [l8ll . 
where the density becomes a highly complicated step func- 
tion under parameter variation, which profoundly simplifies 
the situation. The model was first introduced in 1 38], where its 
parameter dependent diffusion coefficient D{h) was obtained 
numerically, while in |ll|,|39t] the diffusion coefficient for a spe- 
cial single parameter value was calculated analytically. Exact 
analytical solutions for D{h) for all ft, > of this and related 
models were recently obtained in |21]. Since there is a peri- 
odicity with integer values of h, here we restrict ourselves to 
the parameter regime of h E [0, 1] without loss of general- 
ity. In lHHIll] it was found that D{h) displays both fractal 
and linear behaviour, see Fig. \Vlh). To our knowledge, this 
is one of the simplest models that exhibits a fractal diffusion 
coefficient. Being nevertheless amenable to rigorous analy- 
sis, it thus forms a convenient starting point to learn about the 
power of different approximation methods for understanding 
complicated diffusion coefficients. 



HI. CORRELATED RANDOM WALK 

Our first approximation method starts with the diffusion 
coefficient expressed in terms of the velocity autocorrelation 
function of the map, called the Taylor-Green-Kubo formula, 
see |lll|22] for derivations, 

D{h) = lim y / vo{x)vk{x)p* {x)dx 



Kk=0 



VQ{x)p*{x)dx 



(4) 



where p* {x) is the invariant density of the map Eq. ^ modulo 
1, this being equal to one throughout the parameter range as 
we have a family of doubling maps. The velocity function 
Vk (x) calculates the integer displacement of a point at the fc*'' 
iteration. 



Vk{x) = L^^fe+iJ - L^^feJ 



(5) 



In order to create an n*'' order approximation we simply trun- 
cate Eq. (IDl at a given n KM . Hence we obtain the finite sum 



Dn{h) = 



k=0 



VQ{x)vk{x)dx - 



1 



v'^{x)dx , (6) 



which can physically be understood as a time dependent diffu- 
sion coefficient. Looking at how the sequence of _D„(/i) con- 
verges towards D{h) thus corresponds to incorporating more 
and more memory in the decay of the velocity autocorrelation 
function and checking how this decay varies as a function of 
h for a given n. Note that the functional form of Dn{h) for 
finite n is to some extent already determined by our choice of 
integer displacements in Eq. (|5]), however, we have checked 
that for the given model the deviations between using integer 
and non-integer displacements for finite time are minor Sec- 
ondly, we remark that by using this straightforward truncation 
scheme we have neglected further cross-correlation terms that 
do not grow linearly in n, cf. ^. Still, by definition we have 
Dn{h) — > D{h) (n — ^ cxd). Going to lowest order, for n = 
we immediately see that 



Do{h) 



h 
2 



(7) 



which is the simple uncorrected random walk solution for the 
diffusion coefficient ]i2l|]. In Fig. (|2]l one can see that Do{h) 
is asymptotically exact for h ^>- 0. 

Of more interest however are the higher values of n cap- 
turing the higher order correlations that come into play. To 
evaluate these we define a jump function JJ^{x) : [0, 1] — ?> R, 



Jhix) 



fe=0 



(8) 



which gives the integer displacement of a point x after n iter- 
ations. Equation ([8]) can be written recursively as lulll 



ji:ix)^voix) 



J71 — 1 

'h 



jrMM,(x; 



(9) 




D(h) 




Figure 1 : The lifted Bernoulli shift map. A section of the map Mh (a:), Eqs. l|2j and l[3j, is illustrated in (a) for the value of the control parameter 
h — 0.5. The corresponding parameter dependent diffusion coefficient D{h), exactly calculated in I2lll . is shown in (b). 



where Mh{x) is Eq. dU taken modulo 1. This recursive for- 
mula will help when we solve the integral in Eq. (|6]l. Let 

r^(a;) : [0, 1] ^ Rbe defined as 

TJiix) = r r{y)dy, T-\x) := 0. (10) 

Jo 

Using Eq. (|9]) we can solve Eq. JTOl l recursively as 

TJlix) ^ su{x) + ]^T^-^ [Uhix)) (11) 
I 



with 



Sh{x) 



voiy)dy ^ xvo{x) +c, 



(12) 



where the constants of integration c can be evaluated using the 
continuity of T^(a;) and the fact that T^(0) = r^(l) = as 
there is no mean drift in this system. We obtain the following 
functional recursion relation for T'^{x): 



Th{x) - 



irr\2x+h) -\Tj:-\h) 

lTj:-\2x + h-l) -\T--\h) 



\Ti:-\2x-h) 



h-1 ' 



lT--\2x- 



1-h) 



x+ ( 

-lrr\h)-x + {!^ 



2 h 



0<a;< i^ 



2 



<X <^ 

i±^ <a:< 1 



(13) 



r 



Using Eq. ( IT3] ) in Eq. (|6]) via Eqs. ^ and ( fTOl i we can evaluate 
our n*'' order approximation as 



Dn{h) 



rp'n—1 



(h). 



(14) 



So we see that the higher order correlations are all captured by 
the cumulative integral functions TJ^(x). In order to evaluate 
Eq. ( fT4b we construct a recursive relation from Eq. ( fT3] ), 



TJlih) 



n _. 



M|W 



fe=0 



1 

k=l 



\h), (15) 



where 



thix) 



-x + 





h~l 

2 
h+1 



0<x < 



l-h 



1-h 



< a; < 



^<x<l 



(16) 



is sii{x) with the — ^r,"(/i) terms removed. In order to sim- 
plify Eq. ( fTSl ) we write it entirely in terms of Eq. ( fT6] l. Let 



Th{n) 



Z^ 2k-^h 
fc=l 



-fc 



(/i). 



We can write Eq. jTH recursively as 



TH{n)^-T--\h) + -T^{n-l). 



(17) 



(18) 



Substituting Eqs. (fTTT i and (fTsT l into ( fTSl ) we obtain Then substituting Eq. ( fTsT i back into Eq. ( fT9] ) 

(19) 



fe=0 



r,"(/i) 



fe=0 



Ik^h {Mtih) 




A4'W) +^r,(n-l)-ir,(n-l) 



(20) 



r 



we arrive at our final expression 



Tj:{h) 



¥^'^{ 



Mj:{h) 



ri-l 
fc=0 



^i/^O 



MHih) . (21) 



It is helpful to rewrite Eq. ( fTSl ) in the form of Eq. ( |2T] ) as 
it allows us to show that under this method, Dn{h) converges 
exactly to D{h) in finite time for particular values of h, see 
Fig. (|2]i for an illustration. This means that for a specific set of 
parameter values, we can fully capture the correlations of the 
map with a finite time approximation. This convergence is de- 
pendent upon the behaviour of the orbit of the point x = h un- 
der the map Mh{x). In particular, if this orbit is pre-periodic, 
and the values of the points in the periodic loop correspond 
to in Eq. ( fTSl l. then the time dependent diffusion coefficient 
Dn{h) will converge to the exact value D{h) on the n*'' step, 
where n is given by the transient length of the orbit of h plus 
one. For example, let h — 2/5, 

M2/5(2/5) = 1/5 
M2/5(l/5) = 4/5 
M2/5(4/5) = 1/5. (22) 

So h — 2/5 is pre-periodic of transient length one. In addition 



th{2/5) - 1/10 
thil/5) = 
ih(4/5) = 0, 



(23) 



thus to 



for n > 1. Hence we see finite 



^2/5 (%/5(2/5)^ 

time convergence to D{h). This finite time convergence at a 
certain set of points is helpful in understanding the structure of 
D{h) as the fractal diffusion coefficient can be seen emerging 
around these points in the same manner as an iterated function 
system like a Koch curve, see Fig. (|2]). 

Being able to analytically expose the fractal structure of pa- 
rameter dependent diffusion coefficients is the main strength 
of this method. In addition, the convergence of the series 
of approximations is very quick due to the finite time con- 
vergence at certain values of h. Moreover, the fact that one 
only needs to directly put in the map dynamics makes it very 
user-friendly. However, due to the recurrence relation that 



this method is based on, applying it analytically is restricted 
to one-dimensional systems or higher dimensional systems 
whose dynamics can be projected down to one-dimensional 
systems, such as baker maps lUUsl lSTllSSl] . In order to answer 
questions about more realistic, physical systems, one would 
need to resort to numerical analysis. By using families of time 
and parameter dependent diffusion coefficients such as de- 
fined by Eq. (|6]l this is, on the other hand, straightforward, as 
has been successfully demonstrated for many different types 
ofsystemsHiHEHHl. 

This approximation method, represented by Eq. (|6]), was 
criticized by Gilbert and Sanders in I2&1 in two ways; First, it 
was stated that 'this ad hoc truncation has no physical mean- 
ing: if < vqVi >7^ 0, it is not true that higher-order cor- 
relations < VQVk > vanish'. However, there is no assump- 
tion in Eq. (|6]l that higher-order correlations disappear On 
the contrary, this expansion is to be truncated at different time 
steps for exploring the impact of higher-order correlations on 
the convergence of the series by systematically incorporating 
them step by step. Interestingly, as we have shown above, 
there do exist parameter values for this model at which all 
higher order correlations disappear This set, whose number 
of elements becomes infinite for ti — s- cx), holds the key to un- 
derstanding the emergence of the fractal structure in the diffu- 
sion coefficient. There is a clear physical interpretation of this 
set of parameter values in terms of the orbits of the associated 
critical points of the map, as exemplified above. Under pa- 
rameter variation these orbits generate complicated sequences 
of forward and backward scattering, which characterize the 
diffusive dynamics by physically explaining the origin of the 
fractal structure in terms of the topological instability of the 
associated microscopic scattering processes. This physical in- 
terpretation, called 'turnstile dynamics', has been explained in 
detail in iIImIM mi]. 



Secondly, by applying a higher-dimensional equivalent of 
Eq. ^ to the billiard models discussed in ll28ll Gilbert and 
Sanders claimed that in [29*] 'the stationary distribution was 
erroneously assumed to be uniform' . We first clarify that there 
is no room for 'assuming' any stationary distribution in this 
equation. The mathematically exact derivation of the Taylor- 
Green-Kubo formula Eq. (01) is based on time translational in- 
variance of the dynamics, cf. 1 1, 29], and can only be carried 
out if the density p* [x) in Eq. (|6]l is strictly the invariant one. 
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Figure 2: (Colour online) Correlated random walk . In this figure the first four approximations to the parameter dependent diffusion coefficient 
D{h) are illustrated in bold (red) along with the actual diffusion coefficient. In (a) the zeroth order is shown, which is simply the random walk 
solution, in (b),(c) and (d) the first, second, and third order approximations, respectively. At each stage one obtains a set of extrema with linear 
interpolation, which converge quickly to the exact diffusion coefficient D{h). The amount of extrema increases exponentially with n, hence 
we see the fractal structure emerging. 



Hence, there is no choice, and for our map as well as for the 
Hamiltonian particle billiards studied in |28] this density has 
to be the invariant one, which in turn for all these systems 
is uniform. Gilbert and Sanders claim to 'correct this mis- 
take' by deriving a second-order approximation of their bil- 
liard models which is different from the one obtained from the 
method outlined in this section as applied to billiards, com- 
pare Eq. (11) in |281 with Eq. (21) in |29]. Their Eq. (11), 
which they use for their simulations, thus seems to represent 
a mix between the method outlined in this section and the one 
described in the following section. 



IV. PERSISTENT RANDOM WALK 

The next method we look at again starts with the Taylor- 
Green-Kubo formula for diffusion Eq. (HJi. However, rather 
than truncating it, we now approximate the correlations in a 
more self-consistent way by including memory effects persis- 
tently. The key difference to the previous method is that this 
approach models an exponential decay of the velocity autocor- 
relation function beyond the lowest order approximation. This 
method first emerged within stochastic theory as a persistent 
random walk ll30i l3lll and was recently applied to understand 
chaotic diffusion in Hamiltonian particle billiards 128, ,32k, 33l] . 

The main task of evaluating the diffusion coefficient with 
this method is to find an expression for the correlation func- 
tion at the n*'' time step by only including memory effects of a 



given length. We start by defining the velocity autocorrelation 
function as a sum over all possible velocities weighted by the 
corresponding parts of the invariant measure ^* of the system, 

{vo{x)Vn{x)) = ^VQ{x)Vn{x)n* {{vo{x), . . . ,Vn{x)}). 

Vo{x),...,Vn(x) 

(24) 
The different parts of the invariant measure in Eq. ( |24] | are ap- 
proximated by the transition probabilities of the system, de- 
pending on the length of memory considered. These in turn 
are trivially obtained from the invariant probability density 
function p*{x). As a 0*'' order approximation of this method, 
no memory is considered at all, that is, the movement of a par- 
ticle is entirely independent of its preceding behaviour. In this 
case the correlations evaluate simply as 



{vo{x)vn(x)) =0 (n > 0), 



(25) 



thus we need only consider <^Uq(2;)). By Eq. (|4]i the approxi- 
mate diffusion coefficient is obtained as 



Do{h) 



{vl(x))dx 



h 
2 



(26) 



which reproduces again the random walk solution, as ex- 
pected. For the higher order approximations, one must refine 
the level of memory that is used based upon the microscopic 
dynamics of the map. 



A. One step memory approximation 

We now include one step of memory in the system, i.e., we 
assume that the behaviour of a point at the n*'* step is only 
dependent on the (n — 1)*'' step. In ll28ll Eq. (l24l) was eval- 
uated for approximating the diffusion coefficient in a particle 
billiard. For this purpose it was assumed that a point moves to 
a neighbouring lattice point at each iteration. Hence the veloc- 
ity function w„ {x) could only take the values i or —I, where 
I defines the lattice spacing. In order to evaluate the one step 
memory approximation for the map Mh{x), we need to mod- 
ify the method to include the probability that a point stays at 
a lattice point and does not move, hence our velocity function 
can take the values 1, —1 or 0. Let P{b\a) be the conditional 
probability that a point takes the velocity b given that at the 
previous step it had velocity a with a,b G {0, 1,-1}. We 
use these probabilities to obtain a one step memory approx- 
imation. We can write our velocity autocorrelation function 



(wQl'n) = ^ VoVnP{vo)Y\_P{Vi\Vi-l), (27) 

where we let Vk{x) = Vk for brevities sake and p{a) is the 
probability that a point takes the velocity a at the first step. 
We can capture the combinatorics of the sum over all possible 
paths by rewriting Eq. ( |27] | as a matrix equation. 



(wow„) = (01 



Poo 


Poi 


Po-1 \ 


7 


Pio 


Pll 


Pl-1 


p(l) 


P-10 


p-11 


p-1-1 J 


V -p(-i) 

(28) 



where P^a — P(b\a). Equation ( 128] ) can be simplified by 
using the fact that all the paths with a '0' state cancel each 
other out, therefore not contributing to diffusion, and by using 
the symmetries in the system, i.e.. 



P-1-1 = Pll 
P-ii = Pl-1 

p(-l) - p{l) 



Hence Eq. ( |27] l can be simplified to 

{vovn) = ( 1 -1 ) (^ p^"^ p-^ 



-1 



(29) 



p(l), (30) 



which is a simple quadratic form. By diagonalisation the ex- 
pression for the 71*'* velocity autocorrelation function is ob- 
tained to 



(vovn) = 2pil) (Pn - P.^iY 



(31) 



yielding the exponential decay of the velocity autocorrelation 
function < woWn >~ exp(nlog(Pii — Pi-i)) referred to 
above. Substituting Eq. dSTT i into the Taylor-Green-Kubo for- 
mula Eq. ^ by using the fact that p(l) = h/2 gives 



oo _. 

D{h) = Y.{v,v^)- (vl) 



n=0 



hij2{Pii-Pi-i) 



\n=0 



h 
2 



1 - Pll + P] 



(32) 



1-1 



The relevant parameter dependent probabilities can be worked 
out from the invariant density p* (x) of the system and are 



Pii = 





1- 

1 

2 



2h 



0<h< ^ 
l<h<l 



(33) 



and 



0<h< 



Pl-1 



j_ 

2h 



3 

<h<\ 
<h<l 



(34) 



Substituting Eqs. ( [33] ) and ( |34] l into Eq. ( [32] ) we obtain a per- 
sistent one-step memory approximation for the diffusion coef- 
ficient of the map Mh (x). Figure ([3]) shows a plot of the final 
result as a function of the control parameter in comparison to 
the exact diffusion coefficient D[h). 



B. Two step memory approximation 

We now extend the approximation to include two steps of 
memory, i.e., the behaviour of a point at the n*'' step depends 
on what has happened at the [n — 1)*^ and [n — 2)*^ step. 
Let P{c\b, a) be the conditional probability that a point has 
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Figure 3: (Colour online) Persistent random walk approximation . In this figure the first order approximation Eq. l l32t to the exact parameter 
dependent diffusion coefficient D{h) is illustrated in (a), the second order Eq. l |36t is shown in (b). Approximations are in bold (red) along 
with the diffusion coefficient. The major topological changes in the dynamics are picked out by piecewise-differentiable approximations. 



velocity c given that it had velocity b at the previous step and 
a at the step before that with a,b,c (z {0, 1, —1}. For this two 
step approximation, the velocity autocorrelations are given by 



(woWn) 



E 



vav„p{vo,vi)Y\_Pivi\vi-i,Vi_2), (35) 



where p{a, b) is the probability that a point takes velocity a at 
the first step followed by b. Again we proceed by the method 
of 113211 and rewrite Eq. (|35] | as a matrix equation in order to 
capture the combinatorics of the sum, 



(W0«n) 



A" 



(36) 



where r evaluates w„, s evaluates vop{vo , wi ) and A is the 9x9 
probability transition matrix for the system. Unfortunately, 
Eq. ( [36] l cannot be evaluated analytically (see the Appendix), 
so we resort to numerical evaluations. The result is depicted in 
Fig.|3] We see that this method picks out the same topological 
changes in the map dynamics that the previous method did 
and interpolates between them, however, the convergence at 
these points is not as accurate. 

The strength of this method is in modeling the exponential 
decay of correlations that is often found in diffusive systems, 
particularly in Hamiltonian particle billiards [40]. When ap- 
plied to these systems the method is not restricted by dimen- 
sion making it very useful in this setting. However, generating 
by default an exponential decay of correlations is not an ideal 
approach for diffusive systems in which correlations do not 
decay exponentially. In contrast to the correlated random walk 
approach, this method is not designed to reveal possibly frac- 
tal structures of parameter dependent diffusion coefficients. It 
also requires a lot of input about the relevant transition prob- 
abilities, making it unpractical when it comes to analysing 



higher order approximations. In particular, if one was to con- 
sider nonhyperbolic systems [8], or eve n hyperb olic systems 
with less simple invariant measures Ill6l - [l8ll29|] . then deriv- 
ing the transition probabilities of Eqs. ( l33T l and ( |34] | would be 
much more complicated. 



V. APPROXIMATING MARKOV PARTITIONS 

The final method that we look at does not involve the 
Taylor-Green-Kubo formula. Using the framework of the es- 
cape rate theory applied to dynamical systems IIUl2l l34l - [37ll . 
we consider a truncated map Mh{x) defined on [0, L]. By ap- 
plying absorbing boundary conditions to this map, thus gen- 
erating an open system, standard escape rate theory expresses 
the diffusion coefficient in terms of the escape rate from a 
fractal repeller jlSll . Here we use a slight variation of this 
approach by using periodic boundaries. For calculating dif- 
fusion coefficients in simple maps this setting is technically 
easier, because it produces simpler transfer operators than ab- 
sorbing boundaries lllql . We thus consider a closed system, 
whose initial density decays exponentially to an invariant one, 
as quantified by the parameter-dependent decay rate jdec{h)- 
As was shown in II16I - U8I1 . by this modified approach, and in 
complete analogy to ordinary escape rate theory, the diffusion 
coefficient can be obtained to 



D{h) 



hm —^Jde 



■ih). 



(37) 



The decay rate can in turn be calculated exactly if the 
Frobenius-Perron equation can be mapped onto a Markov 



transition matrix. In case of Mh (x) the second largest eigen- 
value xi {h) of this transition matrix determines the decay rate 
according to fSHl] 



ldec{h) = In 



Xiih) 



(38) 



Unfortunately, constructing Markov transition matrices ex- 
actly for even the simplest parameter dependent maps can be 
a very complicated task. Our approximate method starts as 
follows (see lllq, UM for details): For a given value of the 
parameter h, we restrict the dynamics to the unit interval by 
using Eq. ^ modulo 1. We then consider the set of iterates of 
the critical point x = 0.5, which for certain parameter values 
form a set of Markov partition points. This set is then copied 
and lifted back onto the system of size L into each unit inter- 
val. By supplementing this partition with periodic boundary 
conditions, it defines a Markov partition for the whole system 
on [0,L]. The key problem is that the behaviour of the orbit 
of the critical point under parameter variation is very irregu- 
lar Therefore we approximate Markov partitions by truncat- 
ing this orbit for a given parameter value after a certain num- 
ber of iterations. Typically, the resulting set of points will then 
not yield a Markov partition for this parameter value. In order 
to make up for this, we introduce a weighted approximation 
into our transition matrix to account for any non-Markovian 
behaviour. For example, if partition part i gets mapped onto 
a fraction of partition part j then the entry aij in the approx- 
imate transition matrix will be equal to this fraction; see also 
yfl for the basic idea of this approach. 

The motivation behind this method is that at each stage of 
the approximation, whose level is defined by the number of 
iterates of the critical point, there will be certain values of the 
parameter whose Markov partitions are exact. So at least for 
these parameter values we will obtain the precise diffusion 
coefficient D{h), with interpolations between these points as 
defined by the approximate transition matrix. That way, we 
will have full control and understanding over the convergence 
of our approximations. 

We first work out the zeroth order approximation, for which 
we take the unit intervals as partition parts; see Fig.|4]for an il- 
lustration of Mh{x) at system size L = 3. The corresponding 
approximate transition matrix T_{h) is cyclic and reads 
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, (39) 



h 2~2h) 



therefore the eigenvalues can be evaluated analytically lUTl 
Has 

Xi (/i) = 2 - 2/i + 2h cos(2tt/L) 

27r2\ 
2-2h + 2h[l 2" ) {L^oo) . (40) 
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Figure 4: (Colour online) Approximate Markov transition matrix . 
Illustrated here is the map Mh {x) truncated on [0, L] with L = 3 
and periodic boundary conditions. The map is given by the diagonal 
lines (red) and the zeroth order approximation to the Markov parti- 
tion is shown by the thick black lines. The partition parts are simply 
the unit intervals. Note the periodic boundary conditions. The corre- 
sponding transition matrix T{h) is shown below the map. Note that 
this partition is only Markov when /i = or 1. 



as a function of the parameter and length L to 

1 \ h2TT^ 



ldec{h) = In 



1 - h + hcos{2Tr/L) 



L2 



{L -^ cxd). 

(41) 

Using Eq. dTTT l in Eq. dJTl l. the diffusion coefficient is finally 
given by 



D{h) = - , 



(42) 



By combining this result with Eq. ( I38l ). the decay rate is given 



which again yields the familiar random walk approximation. 

The next stage of approximation involves two partition 
parts per unit interval, and for this we simply include the 
critical point x — 0.5 as a partition point. So our parti- 
tion parts are the half-unit intervals on the real line. For the 
next iteration level we include the first iteration of a; = 0.5, 
M^(0.5) — 1 — ft, as a partition point and its mirror image 
about X = 0.5 which is h, and for each higher approximations 
we include one more iterate. However, with these higher ap- 
proximations we no longer obtain a cyclic matrix, so we have 
to resort to numerics to evaluate the decay rate and the diffu- 
sion coefficient. The first three approximations obtained by 
this method are displayed in Fig.|5] 

The main strength of this method is that we know by defini- 
tion where our approximations are going to converge exactly 
in finite time, namely at Markov partition parameter values 
h picked out by each subsequent approximation. In addition, 
the functional form of the interpolation between these points 
highlights areas of self-similarity and therefore gives one evi- 
dence for fractal behaviour even at low-level approximations. 
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Figure 5: (Colour online) Approximating transition matrices . In this figure, the first order approximation to the parameter dependent diffusion 
coefficient DiK) obtained by this method is illustrated in (a), and the second and third orders are illustrated in (b) and (c), respectively, whilst 
a blow up of (c) is shown in (d). The approximations are shown in bold (red) along with the diffusion coefficient diffusion coefficient. We see 
that the functional form of the interpolation in (a) is repeated in (b) at a smaller scale (see the contents of the dashed line box). This functional 
form is again repeated on a still smaller scale in (c) as illustrated in (d). This self-similarity provides evidence that the final function Dih) is 
fractal. 



see Fig. |5] However, this method quickly relies on numeri- 
cal computation and again requires considerable input from 
the user making it unpractical at higher level approximations. 
It seems unlikely that this method can easily be generalized 
to higher-dimensional systems, due to the difficulty to con- 
struct (approximate) Markov partitions and associated transi- 
tion matrices. 



VI. CONCLUSION 

We have studied three different approximation methods 
applied to a particular chaotic dynamical system. For this 
model the exact parameter-dependent diffusion coefficient 
was known beforehand 12111 . Taking it as a reference, the mo- 
tivation was to learn about the capabilities and the weaknesses 
of the individual methods. These are of course not a com- 
prehensive list of the many possible ways to approximate the 
diffusion coefficient of a system, see, for example, 1191 [ill] for 



diffusion in sawtooth and standard maps. However, what they 
do illustrate is the fact that even in our simple one-dimensional 
model studied here, the results that one obtains are very much 
dependent upon the individual method that one uses, and these 
results vary greatly between the different methods. 

By the first method, which relied on a systematic truncation 
of the Taylor-Green-Kubo formula, we saw the fractal struc- 
ture building up fully analytically over a series of correlated 
random walk approximations, as we were able to exactly cap- 
ture the correlations of the system in finite time at certain pa- 
rameter values. This yielded in turn quick convergence to the 
exact results. Using a persistent random walk approach, the 
second method retained an exponential decay of correlations 
even in finite time approximations. However, for the model 
under consideration this approximation yielded convergence 
that was significantly weaker than in case of the other two 
methods. By using a variation of the escape rate approach to 
chaotic diffusion combined with approximate transition ma- 
trices, the third method had our attention focused on areas of 
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self-similarity giving us particularly strong evidence for frac- 
tal structures in the parameter dependent diffusion coefficient. 
This method generated again very quick convergence. Com- 
paring the three different methods with each other demon- 
strates that one is able to tailor the approximate results one 
gets by applying a specific method to the specific questions 
one wishes to answer, or to the specific setting. 

As a side aspect, we addressed recent criticism of the first 
method by Gilbert and Sanders ll28ll . Here we chose a differ- 
ent class of systems than the Hamiltonian particle billiards that 
they considered in their paper. This had the advantage that the 
different approximation methods could be studied more rigor- 
ously. We conclude that the persistent random walk method 
favoured in ll28ll may be more appropriate for dispersing bil- 
liards, because an integral part of this method is modeling an 
exponential decay of correlations, as it is quite common in 
these systems. However, one may question the usefulness of 
this method for diffusive dynamical systems where exponen- 
tial decay is not guaranteed. Here other methods, such as the 
first and the third one discussed above, may yield superior re- 
sults in terms of speed of convergence and identification of 
possible fractal structures in diffusion coefficients. Particu- 
larly the first method has the advantage that it is conceptually 
very simple and quite universally applicable, without making 
any assumptions on the decay of correlations. 



Accordingly, we find the quest for a 'unique' way to ap- 
proximate the diffusion coefficient of a dynamical system, as 
suggested in [28], unnecessarily restrictive. In our view, each 
of the three approximation methods discussed here has, for 
a given model, its own virtue. When one looks to under- 
stand, or display, a particular property of a system and cannot 
achieve this analytically, resorting to one of these approxima- 
tion methods is thus a sensible course of action. 

We finally emphasize that the structure of the diffusion co- 
efficients in more physical systems such as Lorentz gases and 
sawtooth maps are still not fully understood |3, 11]. Particu- 
larly, to which extent these systems' diffusion coefficients are 
fractal remains an open question. Further refining approxi- 
mation methods, such as the ones presented in this paper, to 
highlight areas of self-similarity in parameter dependent dif- 
fusion coefficients in these systems, or to show the emergence 
of fractal structures, would be of great help in answering these 
questions. 



Appendix A: Recurrence relation for the two-step 
approximation 

Equation ( [36] l in a more explicit form is written as 
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which can be simplified using the symmetries in the probabil- 
ities like Pioo = P-ioo and p(l, 0) = p(— 1, 0). In addition, 
if we let TTif be the ij*^ entry of -A" we can use the fact that 



the symmeti'ies of A are the same as A" and reduce Eq. 
to a 3x3 matrix equation. 
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In order to obtain an analytical expression for the matrix in 
Eq. ( lAll i. we would like to obtain a solvable recurrence rela- 



tion, however, this matrix is equal to 
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and unlike for the one-step approximation, a recurrence rela- 
tion is unobtainable, which is due to the introduction of a zero 



state in the velocities. 
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